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Abstract 



■^J- ■ We demonstrate how linear differential operators could be emulated by a quantum processor, 

> 

C**"* ■ should one ever be built, using the Abrams-Lloyd algorithm. Given a linear differential operator of 



m 

00 

to estimate a low order eigenvalue to accuracy 6(l/iV 2 ) is 6((2(5 + 1)(1 + 1/f) + D) log N) 

o 

4h\ 

Q-i, each argument is discretized, v and c are implementation dependent constants of O(l). Optimal 



order 25, acting on functions ip(x\, X2, ■ ■ ■ , xd) with D arguments, the computational cost required 
to estimate a low order eigenvalue to accuracy 0(l/iV 2 ) is 0((2(5 + 1)(1 + 1/f) + D) log N) 
qubits and 0(N 2 ^ S+1 ^ 1+1 ^ U ^ log c N D ) gate operations, where N is the number of points to which 



classical methods require Q(N D ) bits and Q(N D ) gate operations to perform the same eigenvalue 

' 



estimation. The Abrams-Lloyd algorithm thereby leads to exponential reduction in memory and 
polynomial reduction in gate operations, provided the domain has sufficiently large dimension 
D > 2(5 + 1)(1 + In the case of Schrodinger's equation, ground state energy estimation of 

two or more particles can in principle be performed with fewer quantum mechanical gates than 
classical gates. 
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An early motivation for research in quantum information processing has been the sim- 
ulation of quantum mechanical systems The Abrams-Lloyd algorithm 

n n 

instance of quantum mechanical simulation (followed by variations |5||, [y], 
scribe in this paper the application of the Abrams-Lloyd algorithm to estimating low order 
eigenvalues of linear partial differential equations with homogeneous boundary conditions 
(more precisely, Hermitian boundary value problems). The significance of our analysis is two 
fold. First, we generalize the Abrams-Lloyd algorithm to boundary value problems other 
than Schrodinger's equation, which may find application to classical problems. Secondly, 
we quantify computational cost and determine under what conditions we may expect the 
Abrams-Lloyd algorithm to give a reduction in computational work compared to optimal 
classical techniques in order to achieve the same eigenvalue accuracy. 

Very briefly, the Abrams-Lloyd algorithm as originally envisaged for the many-body 
Schrodinger equation is structured as follows. An initial estimate \ip{0)) of the target eigen- 
state is loaded into a multiple qubit register. Controlled application of a unitary operation, 
chosen to correspond to the time evolution operator exp(— iTir) of the many-body Hamil- 
tonian 7i under study for time step r, allows one to generate a sequence of time evolved 
states originating from the initial guess, {^(O)), \ip(r)), \ip(2r)), . . .}. A spectral analysis of 
the sequence of time evolved states recovers the frequency (energy) of the target eigenstate 
(provided the initial guess was "close enough"). The Abrams-Lloyd algorithm is akin to 
a stroboscope for quantum states evolving under a many-body Hamiltonian. If the total 
time of evolution is sufficiently long, while the individual time steps are sufficiently small, 
a high frequency (energy) resolution can be achieved. Following the determination of the 
eigenvalue, the corresponding eigenstate remains in the qubit register. Although the full 
amplitude description of an eigenstate is inaccessible, some information about the state can 
be extracted to a precision ultimately limited by the number of qubits used to represent the 
eigenstate (so for instance, one can test symmetries of the eigenstate). 

The algorithm can be extended to more general partial differential equations rather 
easily. So long as the boundary value problem is Hermitian, we can map our mathe- 
matical problem to a fictional quantum system and apply the algorithm without change. 
The partial differential operator, C, will correspond to a (possibly) fictional Hamilto- 
nian 7i, and an initial guess ip(0) will correspond to an initial wavefunction \ip(0)). In 
other words, quantum mechanical amplitudes represent function values. Less obviously, 
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the sequence of time evolved states, {|-?/>(0)), \ip(r)), \ip(2r)), . . .} has a mathematical ana- 
logue of great use in classical matrix eigenvalue analysis, known as the Krylov subspace: 
spanj^O), exp(— iCr)i/j(0), exp(— i2Cr)ip(0), . . .}. The subspace is generated by repeated 
application of exp(— iCr), although in classical techniques one more typically uses rational 
functions of C. Here, r no longer has the physical meaning of time. Rather, r sets the scale 
for how much phase is applied per application of exp(— iCr). As in the quantum simulation, 
a large total phase applied one small phase step at a time allows a high resolution estimate 
of eigenvalues. We quantify these notions now. 

The computational cost of the Abrams-Lloyd algorithm for a specified eigenvalue accuracy 
is limited as a consequence of three sources of error, expressed here in the language of 
quantum mechanical simulation: 

I truncation error: Discretization is necessary for a computational model based on 
qubits. However, discretization of the continuous problem to N points per coordinate 
results in Q(l/N 2 ) relative error in low order energy eigenvalues due to truncation of 
high spatial frequency contributions. The choice of N must be made appropriate to 
the accuracy that is desired. 

II splitting error: The full many-body evolution exp(— %Ht) over time step r can be 
implemented with universal gates by splitting the full evolution into a sequence of effi- 
ciently implementable unitaries exp(—iH k r), where H = J2 k 'H k - The approximation 
results in an absolute eigenvalue error 0(\\H\\ u 2 +1 t u ), where ||W|| 2 is the maximum 
eigenvalue of the discretized Hamiltonian 7i and v is a constant of 0(1) determined 
by the precise sequence of local operators chosen. Splitting error requires us to use 
small time steps r. 

Ill frequency resolution: A quantum Fourier transform, like any discrete Fourier trans- 
form, can resolve absolute phase to accuracy at best ±n. For a sequence of M 
samples, the relative error in an energy eigenvalue E will be ±it/(MEt). Frequency 
resolution requires us to simulate over a large total time Mr. 
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The optimal way to balance these errors is as follows. Since we are interested in the contin- 
uous problem, we first choose a discretization of iV points per sample so that the discrete 
problem eigenvalue approximates the continuous eigenvalue problem to some desired ac- 
curacy 0(l/iV 2 ). We wish to solve the discretized problem to an accuracy determined by 
the truncation error limit; solving the discrete problem to greater accuracy leads to wasted 
effort since we are interested in the continuous problem, while solving the discrete problem 
to lesser accuracy implies we have wasted effort by choosing too many discrete points N per 
coordinate. We can thereby determine the maximum time step r to keep splitting error no 
greater than truncation error. Next, we can determine the number of time steps M required 
to resolve eigenvalues with the quantum Fourier transform at the truncation error limit. In 
the case of Hermitian boundary value problems, We show in this paper that the resulting 
computational cost is Q(D log N) qubits and 0(N 2 ^ S+1 ^ 1+1 ^ U ^ log c N) gate operations, where 
25" is the differential order of L and c is a constant 0(1). This can be compared with the 
optimal classical cost of Q(N D ) bits and Q(N D ) gate operations. Near optimal classical 
methods approaching these costs do in fact exist |2j|. 

We emphasize that in our analysis, we take a constructive approach wherein we account 
for all the logical operations required to implement the algorithm without recourse to oracles 
that may or may not have physically efficient implementations. This is in contrast to previous 
work including the simulation of spin glass physics ^, and Sturm-Liouville problems (j^ 
and references therin). As stated, our motivation is to compare the computational cost of 
eigenvalue estimation by the Abrams-Lloyd algorithm and optimal classical methods. 

Our paper is organized as follows. In section HI we introduce the one dimensional eigen- 
value problem, which will serve as a useful example with which the principles of the algorithm 
can be illustrated. We derive the truncation error for low order eigenvalues in a way suitable 
for extension to higher dimensional problems. The algorithm itself is described in section 
ITU followed by an analysis of computational cost as it is applied to the one dimensional 
problem in section IIHI A circuit suitable for a 2nd order differential equation is given as a 
concrete example. Generalization of the algorithm to higher dimensional problems is given 
in section HVl along with an analysis of computational cost, where we show a reduction in 
computational work polynomial in N over classical techniques. Concluding remarks about 
the computational efficiency of the Abrams-Lloyd algorithm are given in section |V] 



4 



I. ONE-DIMENSIONAL PROBLEM 



To illustrate the essential features of eigenvalue estimation of differential operators, it's 
instructive to consider a Hermitian one-dimensional problem, which we introduce here in 
some detail. The primary result of this section is a derivation of truncation error in low 
order eigenvalues as a function of discretization. Much of the notation used throughout this 
paper are defined in this section. We begin with a linear, 2S'-order differential operator T> 
that maps a complex valued function ip{x), x [0? 1] to a new function according to the 
rule, 

s 



s=0 v 7 



a (x)il)(x) + J- \ ai(x) ) + • • • 

d s ( d s i/j(x) 



dx s \ S dx s / ' ^ ^ 

where we assume ip(x) has finite derivatives up to order 2S. The coefficients a s (x), s = 
0, 1, 2, . . . , S are finite, real valued functions on the domain x G [0, 1] with finite derivatives 
to order s and satisfy periodic boundary conditions, 

f> = £(l) t = 0,l,..., S (2) 

The minimal smoothness assumed of ao(x) is continuity on x £ [0, 1]. For concreteness, we 
impose periodic boundary conditions upon if)(x) itself, 

^(O) = 0(l) * = 0,1,...,2S (3) 

although more general homogeneous boundary conditions could be insisted upon. Given 
the above definitions, a set of eigenfunctions 4>f{x) with corresponding real eigenvalues A/ 
is defined through, 

2ty/(aO=A/0/(a;), (4) 

and we order the eigenvalues A/, / = 1,2,3, .. . in ascending order Ai < A2 < A2 < . . .. 
The definition and boundary conditions in Eqs. CMH1 guarantee a Hermitian T>, meaning 
Jq dx((j)fV(j)fi — 4>fiV(f)f) = for any pair of eigenfunctions <j>f,(j>f>- All the usual eigen- 
value/eigenfunction properties of Hermitian operators follow. Our task is to estimate a low 
order (/ = 0(1)) eigenvalue A/. 



The most useful expression of the eigenvalue is the Rayleigh quotient, 

Xf — / dx(p*j(x)V(f)f(x) = / dx<p*f(x)C4>f(x) (5) 
Jo Jo 

where we impose unity L 2 norm on the eigenfunctions, 



i \ 1/2 

Li = 



dx(f)*j(x)(f)f(x) J (6) 

in anticipation of the quantum algorithm and the operator C, derived from T> by simple inte- 
gration by parts, is a more convenient (bilinear) operator to work with due to its symmetric 



form, 
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.™*H£^M^ (T) 



s=0 



for any two functions <f(x) and i/j(x). 

It is useful to work not only in the "space" domain x G [0, 1], but also in the "reciprocal 
space" domain of integers, fceZ. The connection between the two representations is defined 
by the Fourier transforms, 



fi 



ipk — / dxexp (— 27rikx)ip(x), 
Jo 

CO 

ip{x) = exp (2irikx)ipk, (8) 

fc=— oo 

where tilde will indicate a reciprocal space representation throughout the paper. Our eigen- 
value Eq. |3] is Fourier transformed to 

oo 

£k,k'(fif,k' = ^f4>f,k, (9) 



k'=— oo 

where, 



C k ,y = Y J ^ k ) S ^s,k^i2mk') s (10) 

s=0 

is the reciprocal space matrix representation of the operator C (and T>). The Rayleigh 
quotient of Eq. El is Fourier transformed to, 

oo 

-V = ^ ^f^k^'^f^', (11) 
k,k'=— oo 

where we now have the Euclidean normalization, 



1/2 



<j) f 



Vfc=-c 



consistent with ||0/||l 2 = 1 and our Fourier transform definition. 

In a classical digital computer, discretization of the domain x G [0, 1] is required so that 
values x can be represented with a finite number of bits. For the quantum algorithm we'll 
be discussing, discretization of the domain will also be required so that values x can be 
identified with a finite number of qubits. We can then sample the spatial domain at the 
points x = 0, 1/N, 2/N, . . . , (N — 1)/N, where N = 2™ requires n qubits. It will be more 
convenient to work with the integers x = Nx = 0,1,2, N — 1. The discrete spatial 
domain of N points allows us to approximate a function ip(x) by a vector, 

^ = (4 N \4 N \---^ N N \) as) 

for computational purposes, where we shall impose Euclidean norm ||V ,( ' iV ' ) Ih = 1- in par- 
ticular, we wish to generate discretized approximations that approach the continuous 
problem eigenvector <pf{x) such that taking an ever greater number of discretization points 
iV gives us the limit linijv-^oo V^^/P = <f>f(x), the factor y/N accounting for Euclidean nor- 
malization of the vector (j)^ and L 2 normalization of the function 4>f(x). We discuss how 
we generate 0^ and how quantify the quality of our discrete approximations as a function 
of N further below. 

In addition to having a discrete approximation to functions ip(x), we shall require discrete 
approximations of the differential operator £ of Eq. [7| in the form of an iV x JV matrix 
acting on vectors Hence, we'll need N x N finite difference matrices, which we 

shall denote to approximate derivatives d s /dx s . There is freedom in choosing finite 

differences to approximate derivatives, here we (arbitrarily) choose the forward difference 
for a concrete example, 

and higher order finite difference, canta generated by A« = (a™)' for integer 
From the very definition of derivatives, we have linijv^oo VN(A^tjj( N ^)x = d s tp(x)/dx s 
if limjv^oo y/Nijj^ = ip(x), the factor \/N again accounting for Euclidean normalization 
of the vector and L 2 normalization of the function ip(x). The subscript arithmetic 
x ± 1 in the definition of finite differences is to be performed modulo- N, consistent with the 
boundary conditions of Eqs. 12131 The resulting matrix operator is, 

^ = E( A w) T ' Difl 8( a 5 JV> )- A w- ( 15 ) 

s=0 
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where (-) T indicates matrix transpose and Diag(-) indicates a diagonal matrix with the 
vector argument along the diagonal. With the above construction for & N \ we can now pose 
a Hermitian matrix eigenvalue problem, 

N-1 

£4>S = AfVS> (16) 

x'=0 

whose solutions will have the desired properties limjv-*oo VN<$p = 4>f{x) and 
lim_A/-->oo —A/, with the obvious restriction / < N. A reciprocal space description 
is useful, for which we introduce the discrete Fourier transforms, 



x=0 
N/2-1 

V 7V k=-N/2 

where u; = exp(27ri/iV) and reciprocal space has been truncated to the set A/" = {k £ Z : 
—N/2 < k < N/2 + 1}. The eigenvectors are assigned unit Euclidean norm in both x and 
/c space representations, 



1/2 



'N-1 

£■ 

E^f»«V /2 = i, (is) 



so that the discrete analogs of Eqs. El IHl are 



JV-l 



*<»> = 

x,x'=0 



r . r, t f,X 

= £ (19) 

which we shall find useful below. 

We shall call — A/| the truncation error, alluding to the fact that we wish to ap- 

proximate Xf with while truncating reciprocal space from all integers Z to the subset 
M . We now proceed to show the well known fact that replacing derivatives by finite differ- 
ences ultimately limits the convergence of A^ to Xf as the number of sampling points N 
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increases. Straightforward application of previously stated definitions gives, 

A^V (iV) J = N (exp (2mk/N) - 1) 

= ^k^(i + e(^\) (20) 



N 2 ) 

where we have made use of series expansions and the fact that \k\ < N/2 to arrive at the 
contribution Q(k 2 /N 2 ). The result holds for higher order derivatives. 

An important parameter in characterizing truncation error is a reciprocal space cut-off 
k((f>f), which can be defined for every <ftf. There will always exist a number k(<f>f) such that 
|<A/,fc| 2 — 0(k~ ( - 4:S+1+ ^) for all \k\ > k(<pf) and some infinitesimal e. This follows simply 
because 0/ must be differentiable up to order 2S, and therefore the series ^fc(2 7r ^) 4S '| ( / , ./',fc| 2 
giving the norm of the 2S th derivative of 0/ must converge. The eigenvalue spectrum of 
the continuous domain operator C is unbounded, and it can be shown that sup{A;(</>/)} does 
not exist. However, since we restrict ourselves to / = 0(1), we can specify a finite k((j)f) 
independent of N. For N/2 > k(<pf), a reciprocal space cut-off k(<j>^ ■) must also exist since 
limjv->oo 0/[fc = 4>f,k- From here on, we shall not distinguish between k(4>^) and k(<pf) as 
the precise value of the reciprocal space cut-off is not needed, but simply its existence. We 
thus define another subset of reciprocal space A4 = {k 6 Z : \k\ < k(<j)f)}. 

We have now collected enough ingredients to find the truncation error \\W-\ f \. We 
assume that N/2 > k((j)f), so that a "reasonable" representation of <pf can be made on the 
discretized domain. By "reasonable", we mean the eigenvalue A/ can be estimated using 
Eq. [TT] and the truncated reciprocal space Af to give, 

X f = E + O (iV-( 2S+e >) (21) 

k,k'eAf 

where the above result arises from the least convergent (highest order derivative) contribu- 
tion to Xj in the region k, k' ^ Af, 

(2mk) s a s , k - k/ {2-Kik') s 4> f>k 

k,k'£N 

= ^ O (k'' {s+1/2+e) k' ( - s+1/2+e) ) 

k,k'^Af 

= O (N-( 2S+ ^) (22) 

where we have made use of |0/,fe| = 0(k~( 2S+1 / 2+ ^) for k > k((f)f). Thus, for N/2 > &(</>/), 
truncation of the reciprocal space sum in Eq. ^2 gives 0(N~^ 2S+e ^) error. 



Using the finite difference error of Eq. [201 the reciprocal space matrix D- N > can be written, 



3 N ] = V 



N 2 



(2irik) s a S:k - k/ (2irik') s (1 + 6 

, .i) 1 

(AT) 



k 2 + k' 2 



(23) 



where we have used the fact that there is some freedom in approximating a s (x) by a s 
We choose to match spectral components af£ = a s ,k, and accept that a££ may exhibit 
oscillation artifacts (Gibb's phenomenon) due to discarding the contributions a s ,k for k G 
Z — M . Note that the smoothness of a s (x), meaning continuity and finite s order derivatives 
for x G [0,1], implies the existence of reciprocal space cut-offs k(a s ). We use Eq. IT^l to 
decompose, 



A, 



(N) _ V" Z*W r (^ + ft ( k ^ 

f ~ 1^ L k',k(p fik 1 + U I — T 



k,k'eM 



+ £ ?«?^(i + e(T)) < 24 > 

k,k'e{Af-M) 



V iv 2 
+_* 

N 2 



where the error summed over M — M. is of order, 

fc 2 + A;' 2 

k,k'e{M-M) 



jy2 / V ^'5+l/2+e^,S+l/2H 



The diminishing contribution of the region M - M to A^ ensures that the relative finite 
difference error Q(k 2 /N 2 ) does not approach unity but remains 0(l/iV 2 ). Collecting the 
results of Eqs. [22 [2H [23 we can express A ^ — A/ as, 

- A / = E *fJPz*j$2 ( 1 + e (^))- E ?/,*^?/.* (26) 

k,k'eN ^ V / / k,k'eAf 

where we have dropped the dependence upon k((f)f) as it shall be of no further use. We note 
that —£ = 0(l/iV 2 ) in the reciprocal space .M, so we can consider 1/iV 2 a parameter of 
expansion in perturbation theory. The lowest order perturbation gives 11 ^ 111 n . • ' 



/ 112 - Wf 

= 0(l/iV 2 ) for a non-degenerate 0/. Degenerate eigenvectors might be perturbed 
substantially, but this is merely the result of there being no preferred basis for the span 
of the degenerate eigenvectors. The same bounds on truncation error can be shown to 
apply to the degenerate case. Noting that Eq. [2S1 is second order in eigenvector and 5(f) / 
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is orthogonal to <pf, the contribution of 5<f)j to the eigenvalue error is 0(1/N A ) and can 
therefore be ignored. The relative truncation error is, 



which is the final result of this section. We emphasize that truncation error arises solely 
from the uniform discretization of the domain x G [0, 1]. 

II. QUANTUM ALGORITHM - ONE DIMENSION 

We present now the quantum algorithm as it applies to the Hermitian, one-dimensional 
boundary value problem discussed in the previous section. We will show the various com- 
putational steps, and the rationale behind them. 

First, we set forth some preliminaries. We will represent a vector with a quantum 
state composed of n — log 2 N qubits whose probability amplitudes are encoded as follows, 



where \x) is an n qubit state storing the binary representation of x. Similarly, the finite 
difference matrix is mapped to an operator, 

N-l 

a^= i*>45<n (29) 

x,x'=0 

and we define the unitary exponential, 



where r is a dimensionless constant whose value is chosen in advance of the simulation and 
where the unitarity of U follows from the Hermitian nature of . The constant r must be 
carefully chosen to arrive at a desired accuracy in eigenvalue without an unnecessarily 
large number of operations. The prescription for choosing r is described further below in 
section IIHI Note that r is now an abstract scaling parameter rather than the time step of 
a quantum simulation. 

We shall call the register of n = log 2 iV qubits the accumulator register. In addition, 
a register of m — log 2 M qubits will be required to count phase steps, which we shall call 




(27) 




(28) 




(30) 
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the index register. Several ancilla qubits will be required, their number depending on the 
desired precision for the coefficients that specify A. The first steps are to load an initial 
state ip( N ^) into the accumulator and to form an equal superposition of all index qubit states, 
giving a complete state, 

1 M-l 

i*> = 7^Ei^ (a0 >i^ ( 31 ) 

V 3=0 

The state \ip( N,) ) is an initial estimate of the field eigenvector of interest. To determine the 
required computational work to arrive at a suitable initial estimate \ip^), it is useful to 
decompose the accumulator state in terms of the initially unknown eigenstates \<fi^), 



f 

\ i — i v-i .v-j 

M 



i*> = ^E E°/i^)w = £«/!*/>' ( 32 ) 

j=0 x =0 f=0 

where acf = ((fr^\i/)^). As will be shown, the probability the Abrams-Lloyd algorithm 
will give an estimate of eigenvalue A/ in a single iteration is |«/| 2 . To obtain an estimate 
of Xf with probability approaching unity, approximately 1/|«/| 2 iterations will be required. 
It is thus necessary for l^^)) to have a large overlap with \(f)f) in order to avoid numer- 
ous iterations of the algorithm. The best technique proposed thus far is that of Jaksch and 
Papageorgiou Q|, where a more coarsely defined 0^°^ is determined first (ie. N < N). Ac- 
cording to the analysis of the previous section, a coarse approximation limited by truncation 
error will allow one to achieve, 



\ a f 



2 < 



2 



5<p 



(No) 
f 



2 



2 



0{l/N%) (33) 



for / = 0(1). Thus, one might solve for a desired <j)f classically (with cost that we will 
discuss later), and load the state 10^°)) into the accumulator with Q(N ) operations. 

We shall now follow the linear portion of the algorithm as it operates on a particular 
component |^/), reintroducing the full superposition over all / in the final (nonlinear) mea- 
surement step. The next stage of the algorithm is to apply the unitary U to the accumulator 
conditional upon the index to produce the superposition, 

1 M-l 

3=0 

M-l 

— ^ex P (yAfV)|0f)|i). (34) 



3=0 
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Only M conditional applications of U are in fact required to form |\E^) from |^/). One 
applies U conditional on j > 1, then one applies U conditional on j > 2 and so forth until 
the (M — l) th conditional U is applied for j = M — 1. The conditional applications of U can 
be performed with a single additional ancilla qubit as follows. With at most logM logical 
operations, one can entangle the index register with an ancilla to form the state \ j)\Cjjr) 
where the ancilla Cjj> = 1 for j > j' and Cyy = otherwise. The j' th application of U 
can be implemented as a U conditional on the ancilla Cjj>. The ancilla is then disentangled 
from the quantum register by running the initial entangling operation once again. 

The operator U acts in the full n qubit Hilbert space of which will in general be 

prohibitively large, but it is nevertheless possible to efficiently generate an approximation 
to U using operations in a few qubit Hilbert space. The structure of is a band diagonal 
matrix resulting from local operations, and thus it has a block diagonal representation in 
the qubit basis of the accumulator. 

To illustrate explicitly some of the key features of the algorithm at work, it's useful to 
consider the simple example where V = d/dx{a(x)(d/ dx)}. The following decomposition 
is appropriate, 



\ 



do 



da 



V 



/ ai 
di 



- N 2 



J 



\ 



a 3 
a 3 



V 



- N 2 



J 







a 2 
a 2 



a 4 

04 



\ a 



A (N,2) A (N,3) 

where d^ = c% + a-x+i- The operators A^ N ' P ^ can be written more compactly, 

A (N,1) = N 2J2^\x)(x\ 

x 

A (N,2) = _ N 2 a w+1 {\x){x+l\ + \x + l){x\} 

x even 

A W3) = _ N 2 a w+1 {\x)(x + l\ + \x + l)(x\}, 



x odd 



where A (JV,1) is diagonal, and AW 2 ) )A (w,3) act 

in one qubit subspaces (conditional upon 

[x/2\ ) in lieu of the full Hilbert space of A^. 
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The unitary U can be approximated to take advantage of the above decomposition in 
several ways. For a general decomposition, 

A^=f A<^, (36) 
P =i 

where for our simple example R = 3, the Baker- Campbell-Hausdorff formulae can be used 
to show, 

R 1 
Un = ] Jcxp (iA (iY ' p V/2) Y\ ex P 

(iA^V/2) 

p=l p=R 

= exp [#V - 1 £ [A^\ [A<"">, A (JV ' fl) ]] r 3 + O (||AW||,r 4 )') (37) 

where terms bilinear in A( n,p >t are suppressed by the symmetry of the product formula 
shown. One may approximate U by Un to take advantage of the efficient implementation of 
exp(iA^ 7V,p V) at the cost of introducing error. The quantum circuit for implementing Un for 
our simple example T> = d/dx{a(x)(d/ dx)} is shown in Fig. ^ The reason for the ease of 
implementing Un is apparent in Fig. [H one applies single qubit unitaries conditional upon 
the evaluation of c%. The Solovay-Kitaev theorem guarantees that the single qubit unitaries 
can be implemented to an accuracy 0(1/7V 2 ) with @(log c A^) universal quantum gates 
We also assume that evaluation of <% requires O(logA^) operations. Roughly speaking, the 
differentiability of a(x) rules out pathological functions that have greater complexity. 

Approximating U by Un implies that the algorithm will give an estimate of the eigenvalue 
of the operator, 

A(f) = AW + o(||AW||'r 2 ) (38) 

instead of the desired eigenvalue \j. We call the error introduced by using Un the splitting 
error, which has value O ^||A^ ||g r 2 j provided ||A^||2T < 1. The splitting error will be 
shown in the next section to limit the computational efficiency of estimating eigenvalues. 
Applying Un rather than U, Eq. (pMjl takes the form 

M-l 

l*/.n> = 7mT. «p(v'4Sr)|^>|j>. (39) 

The eigenvalue is encoded in the phase periodicity of |0/,n)|j), and can be determined to at 
most the ±7r/M precision allowed by a log 2 M bit representation of a full 2tc radians. We 
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FIG. 1: The quantum circuits for applying: (a) exp(iA( Ar ' 1 V/2), (b) exp(zA ( - Ar ' 2 V), and (c) 
exp(iA^ Ar,3 V/2), to the accumulator qubits \x) = IX3X2X1X0) for the decomposition of Eq. (TTH|) 
with N = 2 4 . The ancillae initialized to state |0) are used to store the coefficients d x or a^+i, to 
three bit precision with resolution 5: dx = (d^^ 2 + i2 + d Si o2°) x 5 with d^j G {0, 1}, 
and a similar binary description for a^+i- The circuit "c^" maps |x)|anc) — ► \x)\dx © anc), 
while "a s+ i (even/odd)" maps |x)|anc) — > |x)|a^+i © anc) for even/odd x. Single qubit rotations 
R' p = |0> <0| + exp(i(A^ 2 5r/2)2 p )|l)(l| give the desired diagonal phase shifts for exp(iA^ N '^r/2). 
Single qubit rotations X' p = exp(i(— N 2 5t /2)2 P a x ), where a x = |0)(1| + |1)(0|, to implement the 
desired off diagonal couplings of exp(iA( Ar,2 - ) r/2) and exp(iA^ Ar ' 3 V/2). The values of r and 5 can 
be inferred from the restriction that the operator splitting error is comparable to truncation er- 
ror, described in section ITTT1 The parity shift operators, defined D^\x) = \x db 1), are required to 
shift the block structure of A^' 3 ) so that only operations on the least significant qubit xo need be 
performed. The can be implemented using quantum Fourier transforms (at cost of 0(log 2 N) 
operations) and single qubit rotations. Final disentanglement of ancillae is achieved by a second 
application of a dx" or "<%+i". 15 



briefly review the procedure for retrieving the phase 
the quantum Fourier transform, 



beginning with the application of 



M-l Af-1 

QFT = -= Yl exp (-^ ilm l M ) 10 H, ( 40 ) 
V i=0 m=0 



to the index qubits. The resulting state is, 

M-l 



QFT|^ n ) = E^/lC)IO- (41) 



with the coefficients, 



j=0 

which have square modulus, 



1 M-l 

bl >f = M E ex P (v ( A /> " 2 ^I M ) ) > ( 42 ) 



|2 
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(2tt//M - Aggr) 
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M 2 sin 2 


(2tt//M - \f^T 


)/2 



(43) 



A projective measurement of the index produces \V) where \X^t/2tt — l'/M\ < 1/2M with 
a probability \bvj\ 2 > (M 2 sm 2 (ir/2M))~ 1 . All eigenvalues will satisfy |A^| < tt/t since 
we will impose ||AW|| 2 t < 1 (to be made precise in the next section), so identification of V 
will determine an eigenvalue A^j* uniquely to a precision ±tt/Mt. 

Since we began not with the desired state alone, but with a superposition \^/) = 
Y^k=o a f\^ f)i measurement of the index will determine a particular A^ with relative prob- 
ability |«/| 2 . It is the initial trial wavefunction \ip^ NoS) ) that determines the probability |«/| 2 
of the eigenvalue/eigenvector pair being selected by a projective measurement. 

Upon completion of the eigenvalue readout (via index bits I'), the accumulator is left in 
the eigenstate \(/)f). This is useful since it allows further information to be extracted. For 
instance, one can efficiently test whether \(f>^) has a particular symmetry, such as inversion 
symmetry about a particular point x in the domain. This can serve as a partial check as to 
whether the desired 10^) was indeed selected by the projective measurement. 



III. COMPUTATIONAL COST - ONE DIMENSION 

We now analyze the computational cost for implementing the Abrams-Lloyd algorithm 
for the one dimensional Hermitian problem described in the preceding sections. As pointed 
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out, there are three sources of error that must be considered to determine the required 
number of operations for a given accuracy in eigenvalue estimation. 

First, uniform discretization of the continuous problem to N = 2 n points on the spatial 
domain introduces a truncation error, 



\ f - Af 



A, 



© ( ^ ) (44) 



The truncation error quantifies the accuracy with which the discrete problem represents the 
continuous problem for low order (ie. / = 0(1)) eigenf unctions 0/. To compare algorithms, 
classical or quantum, we may ask how many operations are required to achieve the B(l/iV 2 ) 
accuracy in the solution of the discrete eigenvalue problem. 

Second, splitting into parts so as to approximate U with a product Uu of local 
operators results in what we have termed splitting error. From Eq. EElthe eigenvalue Ay^ 
of Un is, 



,(N) _ ,(JV) 

A /,n - A .f 



+ (||A^>||^ 2 ) (45) 



where we choose r such that ||A^|| 2 t = ||£^||2T < 1. However, from the the finite 
difference formula Eq. ITU and the form of in Eq. [03 the spectral radius ||£^||2 = 
Q(N 2S ). Hence, the splitting error is, 

Agg = A^ + O (iV 6 V) (46) 

which, unlike truncation error, increases polynomially with an increase in the number of 
discretization points N. The splitting error results from the fact that the product Uu creates 
deviations from the true advancement in phase at high spatial frequencies. For example, in 
the system described in Eq. it is the non-commuting nature of advancing even pairings 
of points and odd pairings of points that generates an error with spatial frequency N/2. 

Third, the measurement of phase Ay^r y i a the quantum Fourier transform is limited by 
the uniform discretization of 2n radians into M = 2 m intervals. The limited phase resolution 
allows us to specify Ay^ upon completion of the algorithm to a precision 2tt/Mt. 

The three sources of error allow us to determine the optimal number of index bits m = 
log 2 M, the value of the constant r, and thus the complexity of the algorithm. Obviously, 
there is nothing gained in solving the discretized problem to an accuracy greater than the 
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truncation error 0(l/iV 2 ) if the goal is to study the continuous problem. We can thus allow 
the splitting error 0(iV 6S r 2 ) to be of the same order as the truncation error, 



Since Xfil = 0(1) for our low order eigenvalue with / = 0(1), the phase advancement 
AjJJr < Q(1/N 3S+1 ) for the low order eigenfunction becomes exceedingly small. In order 
to resolve this phase so that our final eigenvalue uncertainty does not exceed the truncation 
error, we require 



thus prescribing the number m = log 2 M of index register qubits. 

The complexity of the eigenvalue estimation can now be stated. The determination of 
a suitable initial guess eigenstate 0j V °' ) requires the determination of an eigenvector of an 
Nq x Nq problem. This can be done classically in Q(Nq) steps, since each of Nq points in the 
spatial domain description of £( N °> must contribute to the eigenvalue. Near optimal classical 
methods are in fact known. In the case of a tridiagonal C^ N °\ bisection gives an eigenvalue 



to 0(l/iVg) precision with 0(iVologiVo) operations 13]. Low order eigenvalues of wider 
bandwidth £( N °) matrices can be determined to the same precision with the same order of 
operations using more complex classical techniques Only a modest Nq is required for the 
probability of a successful iteration of the quantum algorithm, 1 — 0(1/Nq), to be compara- 
ble to unity. Following the construction of an initial eigenstate estimate, this estimate must 
be loaded into the accumulator register, which can be done in Q(N ) steps. We suppose 
that N will exceed Nq by a substantial factor, so that the initial state preparation is a neg- 
ligible cost compared to the remainder of the algorithm. The majority of the computational 
steps in the quantum algorithm are accounted for by the M > O (jV 3 ( 5+1 )) applications 
of Uu, each of which requires 0(log c iV) gate operations for some constant c = 0(1). The 
final quantum Fourier transform requires 0(log 2 M) gate operations, a negligible log(iV) 
contribution compared to the M applications of £/n- 

Thus, to achieve 0(l/iV 2 ) accuracy in the final eigenvalue, at least 0(N 3 ^ S+1 ^ log c N) 
operations and O(-DlogiV) qubits are required. In contrast, an eigenvalue can be found 
using classical techniques to 0(l/iV 2 ) accuracy using Q(NlogN) operations. The quantum 
algorithm requires significantly more work than classical algorithms for the one dimensional 




(47) 




(48) 
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problem. Nonetheless, we show in the next section that the quantum algorithm is easily 
extended to higher dimensional problems where increased efficiency over classical techniques 
is indeed possible. 



IV. HIGHER DIMENSIONAL PROBLEMS 

Here we will generalize the results of the one dimensional problem to the multidimensional 
problem. Many of the arguments presented in the earlier sections are not specific to the 
single dimension domain, and in many cases we can simply replace scalars with vectors. 
The continuous problem we wish to solve involves an operator T> mapping functions 
defined over a .D-dimensional cubic domain x 6 S = [0, 1]® D , to functions 2>0(x). Rather 
than explicitly writing out the general form of a multidimensional Hermitian operator T> 
analagous to the single dimensional operator of Eq. ^ we simply state that T> must satisfy, 



for any eigenf unctions </>/ satisfying T>(j)f = A/0/. We can then define an "equivalent" bilinear 
operator £ that maps any two vector functions ip*(x) and </?(x) to a scalar function ip*C(p. 
This can be done by using the higher dimensional forms of integration by parts, which in 
one dimension allowed us to relate T> to C. We exclude "trivial" problems that are readily 
expressed as a tensor product of single dimensional problems, C = C\ ® £2 <8> • • • <8> C-d- We 
are therefore considering problems whose structure is instead a sum of tensor product terms, 

B 

C = ^£p,l®Cl3,2®...®Cp,D (50) 

13=1 

for some constant B > 1, and where the differential order of each one dimensional Cp^ a is 
2Sp <a . The differential order of C is then 2S = rnaxglJ^ 2Sp j(X \. Of course, we retain the 
Hermitian property 



and the associated eigenvalue/eigenvector properties. Normalizing the eigenfunctions allows 
us to write, 




(49) 




(51) 




(52) 
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which is simply the Rayleigh quotient. 

Discretization proceeds as in the single dimensional case, with each domain coordinate 
%i G [0, 1] discretized to N points. Functions are represented by rank D tensors 
ie. for D = 2 dimensions Vj? is a matrix of numbers, for D = 3 dimensions ^4 is a 
"cube" of numbers and so forth. Partial derivatives are converted to finite differences as in 
the one dimensional case. The operator C is can thus be discretized to a tensor The 
truncation error in the multidimensional problem is, 



which is identical to the one dimensional case because the relative finite difference errors on 
each coordinate are 0(1/A 2 ). 

The implementation of the Abrams-Lloyd algorithm for the multidimensional problem 
proceeds in a completely analagous fashion to the one dimensional case, with the number 
of accumulator qubits Dlog 2 N = Dn so as to represent a volume V = N D . As before, an 
initial estimate of the desired eigenvector </>i is required. A coarse classical simulation can 
produce an eigenvector <pf with Nq < N. Since truncation error scales as 0(l/iV 2 ), the 
required value of N is such that the probability of a successful iteration of the algorithm, 
1 — 0(1/Nq), approaches unity. The computational cost is Q(Nq log N ) classical gate 
operations for generating the initial eigenstate and Q(N®) gate operations to load the state 
into the accumulator. 

The heart of the algorithm is the controlled application of the unitary U = exp(i 
where A^ = \x)C^},(x.'\. As before, U acts within a large Hilbert space, so an approxi- 
mating operator Uu is applied instead. The operator Uu is a sequence of operations acting 
conditionally upon a much smaller Hilbert space than the full D log 2 iV qubits. We quantify 
the size of this Hilbert space now. The multidimensional is no longer represented by a 
band diagonal matrix, but has the structure of a sum of tensor products as in Eq. EOl 

The local nature of A^ is quantified by the maximum number of states |x') for which 
(x'lA^lx) is not zero (maximizing over all possible |x)). This volume, v, is the maximum 




(53) 
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product of matrix bandwidths, 



v = max {(2^,1 + 1)(25 A2 + 1) . . . {2Sp, D + 1)} 




(55) 



where we have used the restriction 2S = maxp^2 a {2Sp a } to arrive at the bound on v. 
It follows that we can split = ^p=i A (7V ' p) where the A (7V ' p) act conditionally upon a 
Hilbert space of r = [log 2 v\ qubits. The size of this reduced Hilbert space is independent of 
domain size N D , so that exp(z t) can be applied to the requisite accuracy (polynomial 
in 1/N) with only G(log c iV) universal gates for some constant c. As before, we assume that 
the function evaluations required for conditional action upon the r-qubit subspace entails at 
most 0(log N) universal gates. The total number of the split up operators A^ N ' P ^ is bounded 
R < v = (1 + 2S/D) D independently of the domain size N D . Thus, Ujj can be applied with 
0(log c N) work for some c = 0(1). 

The approximation Uu can be composed by using a symmetric product as in Eq. |37| so 
that the splitting error is, 



mented |lj,ll5|], uM- F° r the sake of generality, we assume we have a product operator Uu 
correct to order HA^H^t", and set v = 2 to recover the simple symmetric product results. 
In practice, one can not take v arbitrarily large since the number of terms in Uu grows ex- 
ponentially in v. The optimal choice of v is that which minimizes the overall computational 
cost. 

Using the fact that 2S is the differential order of £, the splitting error becomes, 



The final phase measurement through a quantum Fourier transform proceeds as in the 
one dimensional case, with the same precision of ±tt/(Mt) in determining A^, where 
m = log 2 M is the number of index qubits. Requiring that the final eigenvalue be determined 
to the truncation error limit as in the one dimensional case, the same line of reasoning as in 




(56) 




(57) 
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the previous section leads to, 




- y^2(S(l+l/u)+l/u) J 



(58) 



The computational cost of the algorithm is dominated by the M applications of Un, each 
application of Un requiring 0(log c N) number of operations. The computational cost for the 
quantum algorithm is, 



in addition to the cost for finding and loading an eigenstate with coarse discretization N 
along each axis. We assume O(N ). The number of qubits required by the quantum algo- 
rithm is O (log AT). 

We now consider classical costs associated with the multidimensional eigenvalue equation. 
Discretization and reduction of the continuous problem to a matrix equation results in a 
sparse iV D x iV D matrix with a number of bands depending on the spatial derivatives and 
dimensions in the continuous problem. The most efficient and near optimal classical method 
requires 



operations in order to attain a low order eigenvalue with truncation error accuracy 0(l/iV 2 ). 
The method is near optimal in the classical case since the computational cost per each of 
V = N D points in the domain is merely O(logiV). Any classical method must "visit" 
each point in the simulation domain in order for that point to influence the outcome of the 
classical calculation, hence the classical computation cost is Q(N D ). Of course, the number 
of bits required is Q(N D ). 

The maximum improvement in computational efficiency provided by the quantum algo- 
rithm presented is, 



with respect to the best known (near optimal) classical algorithm. From the above, we see 
that the domain dimension must satisfy D > 2{S + 1)(1 + 1/u) in order to see any improve- 
ment using the Abrams-Lloyd algorithm. In particular, we have S — 1 for Schrodinger's 



K Q = 0(M \og c N) = 0(N 2{s+1)il+1/u) \og c N) 



(59) 



K c = 0(N D log N) 



(60) 




(61) 
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equation and we can identify D/3 as the number of particles in space (3 degrees of freedom 
per particle, neglecting spin). A many-body eigenvalue calculation is more efficient than 
classical simulation for particle number D/3 > (4/3) (1 + For the case where Uu is a 

simple symmetric product, v = 2 and we require D/3 > 2 in order to see improved compu- 
tational efficiency Higher order approximations, v > 2 will result in two (spinless) particle 
calculations already being done more efficiently using the Abrams-Lloyd algorithm. 

We now discuss the generality of the results for domains other than the simple hypercube 
S = [0, 1]® D discretized to V = N D points. A more complex domain S' can be had by 
deleting regions from S along planes defined by the uniform discretization scheme. The 
computational cost incurred is that required to ensure the probability amplitudes in \ip) do 
not "spill" into the deleted regions S — S'. This is easily done by circuits such as those in 
Fig. Q wherein quantum gates can be used to determine the conditional application of a 
few-qubit operator through out the simulation domain. The computational cost is therefore 
proportional to the classical cost of determining whether a point x is in or out of the specified 
domain S' subtended by the hypercube S. As an explicit example, the subcircuit <% + i of 
Fig. Qfor applying exp(iA^ Ar ' 2 V) of Eq. EH] can be made to compute |af)|anc) — > |x)|anc) 
for x G S — S' and |x)|anc) — > \x)\a,x+i © anc) for x G S' . The effect of this operation is 
to conditionally apply exp(iA^ Ar,2 V) to those points x G S'. Clearly, S' = S is the simplest 
domain to consider as there is no added computational cost, but more complex domains are 
accessible at only the modest cost of describing the domain with a Boolean function. 



V. CONCLUSION 



Our analysis of the Abrams-Lloyd algorithm raises several questions. Firstly, it is natural 
to ask what sort of qubit phase rotation accuracy is required during the application of Uji to 
the initial guess eigenstate. The phase that is applied to qubits by the operator Uu during 
the computation is of the same order as the phase applied to the highest order eigenvector: 
A^qT where the eigenvalue = <d(N 2S ) for a differential operator of order 2S and 

t = Q(1/N 3S+1 ) for a second order splitting formula. The magnitude of the phase rotations 
applied to qubits is therefore Q(l/N s+1 ). The relative accuracy with which the phase must 
be applied is 0(l/iV 2 ) if the final eigenvalue estimation is to be accurate to the truncation 
error limit of 6(1/A 2 ). Thus the absolute accuracy required from single qubit rotations is 
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fl(l/N s+3 ), independent of the number of dimensions D. The absolute accuracy is a small 
quantity for very modest values of N = 100 (representing a relative eigenvalue accuracy of 
10~ 4 ) with a second order operator (2S = 2). Angular resolution of 10~ 8 in the control of 
qubits represents a technical feat, but thankfully the principles of fault tolerant quantum 
computation jl?], I3| can be applied here to lessen the accuracy requirements for physical 
qubit operations. 

Secondly, it is tempting to compare the quantum and classical algorithms for the simu- 
lation of dynamical evolution. The Abrams-Lloyd algorithm simulates the dynamics of the 
Schrodinger equation dip/dt = Dip for some (possibly fictitious) Hamiltonian represented 
by £>, but only limited detail of the dynamics in a quantum simulation are accessible. The 
probability amplitudes characterizing a register of D log iV + log M qubits can result in at 
most D log iV + log M classical bits of information being extracted by measurement (by the 
Holevo bound). For instance, in order to obtain the eigenvector coefficients 4>^\ at least 
0(iV D /logiV) iterations of the algorithm would be required. This is in contrast to a clas- 
sical simulation of dynamical evolution where 0(A) bits would be required to store a state 
at a single dynamical step, and Q(NM) bits are required to store the entire evolution of 
an initial state over M dynamical steps. We emphasize that the strength of the Abrams- 
Lloyd algorithm is not in its ability to provide great detail into dynamical evolution but 
rather in extracting useful classical information (such as eigenvalues) from a very compact 
representation of that dynamical evolution. 

Finally, the analysis of the Abrams-Lloyd algorithm raises the question as to why the 
eigenvalue convergence for low dimensional problems (ie. small D) is less than that of opti- 
mal classical approaches. Part of the answer lies in the classical theory of matrix eigenvalue 
calculation. An important tool for numerical estimation of eigenvalues is the Krylov sub- 
space, which is defined to be the span of the set {ip, Aip, A 2 ijj, . . . , A M '~ 1 vp} for some initial 
guess vector i/j, some hopefully small constant M' < N D , and some N D x N D matrix A of 
which we seek several low order eigenvalues. The Krylov subspace is spanned by at most 
M' vectors, rather than the full iV D vector space of A, and so projecting onto the Krylov 
subspace gives an efficient means of estimating eigenvalues/eigenvectors of A. If the matrix 
whose lowest eigenvalue is sought is A^, then we might choose A = — fil)^ 1 where \x 
is an initial estimate of the eigenvalue sought (the eigenvalues of being simply related 
to those of A). With A = (A^ - the vector A j ip converges exponentially towards 
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the eigenvector (jy. ' whose eigenvalue minimizes — n\. In contrast, if A = exp(iA^ N ^) 
as in the Abrams-Lloyd algorithm, there is no such convergence towards a target eigenvector 
since the eigenvalues of A are of unit norm. The unitarity of quantum gates restricts eigen- 
values to lie on the unit circle in the complex plane, which is a poor eigenvalue distribution 
from the perspective of estimating a target eigenvalue [19]. This leads to the question of 
whether controlled decoherence can be used to produce non-unitarity evolution to accelerate 
the selection of a target eigenvector with a net reduction in gate operations/delay. 
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